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Abstract 

A new method for interface tracking is presented. The interface representation, based 
on domain decomposition, provides the interface location explicitly, yet is Eulerian. This 
allows for well established finite difference methods on uniform grids to be used for the 
numerics of advecting the interface and other computations. CFL-stable and second order 
accurate explicit and implicit time-stepping methods are derived. Numerical results are 
given to substantiate stated convergence properties, as well as convergence in interface 
curvature and mass conservation. Our method is applied to a boundary integral formula- 
tion for Stokes flow, and the resulting integrals are analyzed and treated numerically to 
second order accuracy. Finally, we embed our method in the familiar immersed boundary- 
and immersed interface methods for two-phase Navier-Stokes flow. 



1 Introduction 

Interface tracking has received considerable attention in the last 20 years. Much of this has 
been motivated by applications in fluid mechanics, e.g. multiphase flows and crystal growth. 
There is also strong interest in the mathematical community on a more general class of 
problems that involve studying PDEs with dynamic free boundaries, so called free boundary 
problems. The present work provides an approach to interface tracking that differs at a basic 
level from the established methods. 

An interface may be thought of as a boundary between two or more sub-domains - in 
represented as one or more curves, and in represented as surfaces. In the language of flow 
applications, we let u{-,t) denote a convective field that deforms the interface according to 
the ODE 

f = u(r,,), (1) 

where T denotes the interface, e.g. given in 2D by a parametric curve, T = {{x{s,t),y{s,t)) : 
t > 0,s e [a b)}. The task of any interface tracking method is to provide a concrete represen- 
tation of r and an equation of time evolution equivalent to (1) that lends itself to practical 
computation. 
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There are several approaches present that provide such a representation, e.g. impHcit 
methods (inchiding level set, phase field and volume of fluid methods) and pure Lagrangian 
methods (such as front tracking) . These are typically introduced in the context of a particular 
application domain, where additional questions often arise. We shall not attempt to survey 
this entire field. 

The level set method, introduced by Osher and Sethian in [1], and with subsequent work in 
e.g. [2, 3, 4, 5], is widely used. This is a fundamentally Eulerian method where the interface 
is not exphcitly specified. Instead it is defined as a contour, or level set, of a function 
which is defined over all of the embedding space. This provides several advantages, including 
that it generalizes well from 2D to 3D, and that topology changes can occur without special 
treatment (though this feature is hard to control since it depends directly on the discretization 
resolution). Notable drawbacks include the fact that the interface tracking problem is now 
discretized in one dimension higher than the interface itself (at a computational cost), and that 
a sharp (explicit) interface representation may be desired in various applications (including 
e.g. surfactant problems [6, 7] and various methods for enforcing interface jump conditions, 
such as the immersed interface method [8, 9, 10]). Other implicit methods, such as phase 
field [11, 12] and volume of fluid (VOF) methods [13, 14] are widely used in computational 
mechanics. 

The front tracking methods introduced by Unverdi and Tryggvason in [15], on the other 
hand, discretize (1) directly as a discrete set of points which are evolved individually by 
integrating the ODE. Subsequent work can be found in e.g. [16, 17]. Front tracking methods 
preserve the dimensionality of the interface tracking problem and provide the interface location 
sharply. An explicit parametrization of the interface may none the less need to be constructed 
(e.g. as in Glimm et. al [17]) to deal with various issues, such as point distribution or 
computation of interface properties, e.g. curvature. 

Recently so called hybrid methods have been introduced that apply ideas from level set, 
front tracking and VOF methods to each other in various permutations. See e.g. [18, 19, 
20] and the references therein. Whereas these methods have successfully remedied some 
flaws of the basic methods, they clearly add practical and theoretical complexity to already 
complicated methods. 

In this paper we take a fundamentally different view on the representation of the interface 
and its dynamics that provides a method that is both Eulerian and explicit. In essence, 
we decompose the interface into a collection of segments, each of which is described as a 
single valued function. Our inspiration for this comes partly from two established fields: 
domain decomposition (see e.g. Demmel [21, Sec. 6.10], and Toselli and Widlund [22]) and 
differentiable manifolds (see e.g. Nakahara [23]). This work is a reformulation and extension 
of the segment projection method by Tornberg and Engquist [24, 25]. 

In addition to the conceptual distinctiveness of this approach, which we think is relevant in 
the somewhat mature field of interface tracking, our method provides some new possibilities 
that will help tackle complex applications. Surfactant problems (i.e. where additional PDEs 
are solved on the interface) is one such example, and contact line problems (i.e. where the 
interface interacts with a solid wall in the presence of some flow) is another. For the latter case 
it becomes possible with our method to pose boimdary conditions on the interface itself in a 
natural way. We also note that this method is efficient and practical to implement, at least 
in 2D. In 3D the efficiency improvement over other methods is noteworthy, but technical 
challenges with representing general closed surfaces have prevented us from deploying our 
method in 3D in full generality. 
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Figure 1: Interface split into segments (ends marked with o), with two segments drawn bold, 
and their local representations. 

In Section 2 our method is introduced and the relevant equations are derived. Then, in 
Section 2.1, we present an explicit and implicit numerical method for our equations, including 
the domain decomposition approach and various related components. We give numerical 
results for the interface tracing method in Section 3 and conclude with two applications: a 
boundary integral method for interface Stokes flow (Section 4), and immersed interface- and 
boundary methods for two-phase Navier-Stokes flow (Section 5). 



2 Interface tracking method in 2D 

Let a simple plane curve be represented by T C M^. Introduce a parametrization, s G [a 6), of 
r at time t, T = T{s, t) = (x(s, t), y{s, t)), as in the previous section. We will assume that the 
interface has sufficient regularity, e.g. with x and y ^ C'^. The evolution of T is described by 
the ODE (1). We will view the the interface as set of segments, ji C T, U'ji = T. Each such 
segment we will choose such that it can be expressed as a single valued function, as presented 
in Figure 1. That is, assume there exists a in interval of s' C s and mapping M. such that 

Mj = {C,f{0) (2) 

with 7 := r(s'). We call this representation, Ai^ := 710c, the segment in local coordinates, 
as opposed to 7, which we will refer to the segment in global coordinates. The mapping will 
be composed of a translation, b, and a unitary rotation, T, of some angle 9. We define this 
mapping and its inverse as 

A^x:=T(e)x + b (3) 
£x:= M-^x = T{-e){x-h), (4) 

for any x E M^. The inverse mapping £ will naturally play the role of transforming the 
segment from local to global coordinates, i.e. 

>C : 710C 7- 



3 



Now we proceed to providing an equation for the time evolution of each segment. Differentiate 
the transformation (3) to obtain 




Hence, we have 




which we may write more conveniently as 



dt Vdi 



0. (5) 



We shall consider first the case where the mapping parameters 9 and b are both constant and 
then, in Section 2.8 we present an extension where we have them as time dependent functions, 
e{t) and b(t). 

For the first case it is evident that 

=r(e)u(7,t), 

where J is the Jacobian and we have used the interface ODE (1). Then we have the segment 
advection PDE given by 

% + v{l.t)^^=w{^,t), (6) 

where we have introduced T{0)n{^, t) =: {y, w). This is an advection-type PDE with a source 
term. Note that the parameters v and w depend on the interface location, so the PDE is 
non-linear. In Section 2.1 we give an explicit and an implicit method to treat this PDE 
numerically in a simple manner. 



2.1 Segment discretization 

We now present numerical methods for treating the segment advection equations (6). There 
are two primary issues to address. First, as previously noted, we deal with non-linear equa- 
tions. However, these are hyperbolic - in quasi-linear form - and the non-linearity presents 
no major concern. 

The second complication concerns boundary conditions and is more interesting: there 
are no boundary conditions at the ends of each segment end if the interface is closed. We 
have decomposed the interface into segments, which leaves us with the requirement that 
segments match where they meet and that smoothness is preserved there. It is then natural 
that we devise a domain decomposition method that solves the segment advection equation 
concurrently in time for all segments. 
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We discretize the segment functions, /(^) on uniform grids in ^, e.g. = a + iA^. 
This constitutes a simphfication over other methods - we can use difference formulas over 
ID uniform grids for our numerical treatment of the segment advection PDE as well as 
all invariant interface quantities that we may need (such as curvature, normals, quadrature 
formulas etc.). 

The convective field, u(-,t), may be an explicit function, as in our test problems later 
in this paper. We also give examples later where u is a functional over the interface (in a 
boundary integral method), and, finally, a numerical solution operator to the incompressible 
Navier Stokes equations. Prom the point of view of the interface tracking method, u(-,f) is a 
black box that orders a convective field to the grid function on ^. 

In anticipation of these more general coefficient operators it is sensible to introduce Strang 
splitting for the pair (/, u): Let u be discretized on staggered time levels, i.e. that 

U"-V2 ^ un+1/2 J jn 

(7) 

It is clear for linear systems that this method maintains second order accuracy, provided 
second order accurate time step schemes for / and u (if u is governed by a differential 
equation) . 



2.2 Lax-WendrofF method for segments 

For the explicit method we will derive a Lax-Wendroff method with central difference approx- 
imations in space (see e.g. the textbook by LeVeque [26]). The same method was used by 
Tornberg et. al. in the segment projection method [24], and we present it here for complete- 
ness. Consider 

ft + v{^{t),t)f^ = w{^{t),t) (8) 

and let At represent the time step. Note that at a particular point in time, v and w are 
functions of ^. We have a Taylor expansion 

fit + At, = fit, + At/t(i, 6 + -^fttit, 6 + • • • 
where the task is to replace time-derivatives with derivatives in space. From (8) we have 

ft = ^^;(x,^) - vix,t)f^. 

This implies that 

fit = {wix,t) - v{x,t)fi:)t 

= ^w;(x, t) - f^^vix, t) - vix, t) (^"^(x, t) - f^^vix, t) + i;(x, t)f^^^ . 
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Plugging into the Taylor expansion gives the semi-discrete scheme 



/(e, t + At) = fiC, t) + At (w;(x, t) - t;(x, t)/^) + 



+ 2 



^x,i)-/^^t;(x,t)). (9) 



At^ ( d 
,di 

Recall that x corresponds to / via the mapping M.. 

For the discretization in space we introduce the usual difference operators Dq, -D- and 
i.e. for some grid function on the ^-grid, discrete derivative approximations are obtained 
via 

(L»og)i = ^(gi+i -gi-i), = -gi-i), = ^(gi+i - gi). 

Despite the fact that this Lax-Wendroff method uses centered difference approximations to 
discretize a convective equation, it is stable provided that the CFL condition is satisfied [26] . 
We note that we will maintain formal second order accuracy even if the time derivatives in the 
last terms are only approximated to first order and, hence, we use a forward approximation 
here by taking an Euler step to get an approximate interface location where the coefficients 
can be reevaluated: 

/* = /" + At(u;(x", r) + ^(x", i")L'oD (10) 

so that 

At )■ ^^^^ 

We emphasize that the evaluation of the coefficients are typically done in global coordinates, 
i.e. via the mapping and its inverse. See Section 2.7 for a summary of the method. As 
previously noted, one may often (depending on how u is computed) want to apply Strang 
splitting (7), i.e. taking all time levels above to t"'^^^'^. 



2.3 How to close the system 

It remains to discuss the how to couple the segments together. With an equation of the 

form (8) it is clear that information will be propagated along characteristics across segment 
boundaries. In the language of numerical methods for conservation laws, any difference ap- 
proximation must then look upwind across to the adjacent segment at any boundary where 
there are incoming characteristics. 

In our proposed method above, centered differences are used. The Lax-Wendroff method 
introduces a diffusive term which stabilizes it (see LeVeque [26, pp. 100-102]). Numerically, we 
then have information propagating in both directions at each segment boundary. In practice 
this means that we need ghost points at each end of each segment. 

To evaluate Dq we can exchange ghost values at time t". Note that adjacent segments 
have different parametrization, so one may need to construct a small interpolant from the 
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Figure 2: Ghost point (x), evaluated from part of adjacent segment (solid line), with local 
segment grid function (dots). 

end of the adjacent segment in order to evaluate the ghost point, see Figure 2. That is, let 
(^fc, fkiO) be the local representation of segment k (or part of). To view this segment in the 
coordinate frame of segment j, simply use the mappings: (^j, fj{0) = ■M.jCkiCk, fk{S,))- With 
this exchange of ghost points, the system is closed in the sense that the concurrent updates 
of all segments is equivalent to (1). 

From a theoretical point of view it is important that there exists a finite overlap between 
segments (for the manifold properties to exist). That is, one must be able to go from one 
segment to another via the global coordinate representation. If this property is not satisfied, it 
becomes impossible to evaluate the ghost points as discussed above. Thus, it will be necessary 
as the interface deforms to expand and/or contract the domain of definition of some, or all, 
segments. However, growing a segment involves an adjacent segment in the same fashion as 
the ghost points. We emphasize that this is handled on a regular grid and, thus, does not 
pose a particular difficulty (as one may view the point distribution problem in front tracking 
methods) . 

2.4 Implicit methods for segment advection 

The natural framework for treating hyperbolic PDEs, including the segment advection equa- 
tion ft + v(S^,t)f^ = w{S,,t), are explicit finite volume methods (to which the Lax-Wendroff 
method may be said to belong). None the less, implicit methods are often proposed for 
interface tracking problems, with enhanced time-stability in mind. 

As a motivating example, consider a multi-phase flow application. A common approach is 
to split the tasks of interface tracking and solving flow equations into interleaving steps. It is 
evident that numerical schemes for both interface tracking and bulk flow should be governed 
by CFL-type time-step restrictions, and indeed, a wealth of numerical methods exist with 
this property. However, it has been observed, e.g. in Peskin and Tu [27], LeVeque and Lee 
[8], Le et. al. [10] and Li and Ito [9, p. 226], that this split leads to a stiff system that is 
no longer stable in the CFL regime. Thus implicit methods are advocated to maintain the 
stability region. 

A time-discrete implicit method may take the form 

which naturally suggests Newton- or fix-point iteration methods to solve for Z""*"^. In our 
framework it is natural to iterate the segments concurrently, in lieu of the coupling between 
adjacent segments (previous section). Let the leading super-index of f'^'"^ refer to segment k 
in the collection that constitutes the interface, and let the latter index denote time-level or 
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iteration count. Then an iterative method is outlined as follows: Until convergence, ||/'^'('+^) — 
/*^'«|| < e,VA;, do 

1. Exchange ghost points between all adjacent segments, 

2. For each segment, k, compute new iterate, e.g. 

for some Q that defines the iterative procedure, 
and let f''''^+^ be the converged iterate for each segment. 

2.4.1 Crank-Nicholson (CN) method for segments 

We now proceed to a concrete method along the lines of the previous section. In [8] and [10] 
a Crank-Nicholson method for (1) is suggested, 

^n+l ^ + ^ (u(x",r) + U(x"+\r+l)) , 

together with an iterative method for solving this nonlinear system. This suggests that if we 
have 

^ = Q(i)/(x,i) + P(t) (12) 
where Q is (possibly) a spatial differential operator, we should study methods of the form 

jn+l ^ jr(/n+l^ .- fn^^ (^QU p ^ Q^+l f n+1 ^ pn ^ _ ^^^^ 

Prom (6) we get a semi-discrete equation on the form (12) with 

Q = -v{-K,t)DQ 
P = w{yL,t). 

We use the following iterative procedure to solve (13) for f"'^^, which follows the approach 
presented in [8]: 

1. Initial guess, /(°) = /" - Atv"Do/" + Atw". Set i = 0. 

2. Evaluate u^*-*. This might involve solving flow equations (or more). 

3. Exchange ghost points between all adjacent segments. 

4. Evaluate the residual, 

5(/«) = /(^^ - r - ^{-v^inD^r - ^«(r+^)z?o/(')+ (14) 

+w''{e) + w^^{e+^)) 

5. If ||5(/^*^)|| < e then terminate and set f^^^ = f^^\ Else increment i, compute new 

and go to step 2. Computing the next iterate should be done cleverly - LeVeque[8] 
suggests a BFGS or SRI method. 



The BFGS method amounts to a Newton-type iteration where, after applying the Sherman- 
Morrison formula, an approximation to the inverse Jacobian (denoted B) is computed con- 
currently with the solution to the system (from [9, Sec. 10.2.5]): 



where 



T 



Here, as in the explicit method, the coefficients f and u)*-'^ are evaluated from /(*^ in global 
coordinates via the mapping (3) and its inverse. We have found that given the initial guess 
computed from the previous time step, this iterative method is reliably convergent in two to 
five iterations, which is well within the practical range, even if u is expensive to evaluate (see 
Section 5.4 for more remarks on this). 

In practice, we have found the iterative method to relax the CFL time step restriction 
present in the explicit method. However, arbitrarily large time steps can not be taken, in part 
due to the requirement that the segments maintain a finite overlap throughout the iterative 
procedure (for the evaluation of ghost points). More remarks on the stability of these methods 
when applied to the Navier-Stokes equations are given in Section 5.4. 



2.5 Error analysis 

The truncation error in each of the methods above is 0{A^'^) + 0{At'^) + 0{NtA^P), where 
p is the order of the interpolation method used to evaluate the ghost points from adjacent 
segments and Nf := T/At is the number of time steps taken. One can show, by constructing a 
partition of unity, that the accuracy order in the interface itself is equal to the order obtained 
for each segment (i.e. second order). 



2.6 Finding and maintaining segments 

An interesting question naturally arises from the assumptions on our method: does there 
exist a decomposition of a curve into segments, and if so, how do we find it and the mapping 
parameters in (3)? 

In a continuous setting it is clear from elementary calculus that this partitioning exists 
provided sufficient regularity on T. However, the curve may have, or develop, a interval 
with arbitrary curvature. This provides a problem in the discrete setting because it implies 
an upper bound on the curvature that can be reasonably resolved with a given segment 
discretization. There is also the question of optimality: partitioning the interface with as few 
segments as possible. We do not address the latter issue, because there is no benefit for our 
method in having very few segments. 

We will also give some remarks on how to verify that the segments, as the interface 
deforms, remain valid (i.e. single valued). But fist the partitioning problem. 
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Figure 3: Winding angle, positive with respect to the orientation indicated. 



2.6.1 Partitioning the interface 

Given a discrete plane curve {xi,yi), i = 1 . . .N, the task is to split the index space into 
intervals that are monotonous with respect to some basis vector. To be precise, consider 
an interval i = p,p + 1, . . .p + M and associate with this segment a mapping (3), i.e. M : 
{x, y) — /). The segment of x from p io p + M \s strictly monotonous with respect to the 
direction defined by if ^j, z = 1, . . . , M, is strictly monotonous. 

A large body of work exists in the field of image processing that aims to reconstruct geo- 
metric features from image data with robust methods. Some related work to our partitioning 
problem can be found in [28, 29, 30, 31]. However, we were not able to find a suitable method 
there or in the references therein. 

The method we propose uses the winding angle, A9i, as defined in Figure 3, and a cumu- 
lative winding angle, 

M 

9i ■.= Y.A9i. 

Then we may say that an interval q = {i = p . . .p + M} is monotonous if 

max 9^ — min 9^ < vr. (15) 

In practice it makes sense to have segments that sweep less than the maximal winding angle, 
TT. So we we iterate over the curve, {xi,yi) from some starting index p and break to a new 
segment when we encounter an M such that max 9^ — min 9^ > r]. Here 77 < vr is a parameter 
that affects how many segments we obtain, e.g. 77 = 7r/2 gives a partitioning of a circle into 
four segments. 

2.6.2 Segment validity 

It is clear that, if the interface deforms sufficiently, the collection of segments initially con- 
structed may fail to correctly represent the interface. That is, given an interface, it may not 
be possible to satisfy the constraint of having single-valued segment functions with a given 
number of segments. Hence, to handle cases with large deformation, we need to monitor the 
segment description and recompute the partition if necessary. We choose to put an upper 
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bound on the variation of the discrete segment function, 

l/n+l — fn\ < K\^n+1 — Cn|) 

where if is a parameter. Also, it is sensible to require that the segments be at least d grid- 
points long. If the interface develops a point of very high curvature, the length condition will 
suggest that a finer grid is needed to resolve the interface. Typical values may be if = 5 and 
d = 20. 

2.7 Summary of method 

Now we summarize the proposed method: 

Time step (explicit method) To advect the interface from tn to tn+i do, 

1. Exchange ghost points between all segments, see Section 2.3. 

2. For each segment (individually): evaluate coefficients ?;(x"', i"), it;(x", t") by using 
the transformation to global coordinates. Then, take Euler step (10) to obtain 

/*, and reevaluate the coefficients v{x.* ,t'^^^),w{-x* ,t'^~^^). Finally compute the 
correction step (11), getting /""^^. See Section 2.2. 

3. Verify that segments are valid, see Section 2.6.2, and possibly grow/contract seg- 
ments so that the evaluation of ghost points possible in the next step. 

4. If segments are not valid, compute a new partitioning (Section 2.6.1). 

The implicit method is proceeds similarly, the key difference being that item two (2.) above 
is replaced by the iterative procedure given in Section 2.4.1. 

2.8 Velocity field decomposition and time-dependent mapping 

Finally, before presenting numerical results and applications, we make an extension of the 
method, where the mapping parameters are allowed to vary over time. Is essence, this poses 
the method hitherto presented in moving reference frames. 

We see in a number of interface tracking applications that there often exists large scale 
deformations that translate the interface or rotate it as a solid body - physical flows are e.g. 
often directional with some length scale. These kinds of motion may be easily treated in our 
method by simply updating the segment mapping parameters (the offset b and the rotation 
6), whereas interface deformations are treated by solving the advection PDE (5), see Section 
3 for numerical results. To make use of this, consider a decomposition of the velocity field u, 

U(X, t) = Utrans -|- cD X (x - b) -|- Urest(x, t), (16) 

where cD = ljz. Since we know Uj = {vi,Wi) at many points Xj we may solve for the decom- 
position parameters co and Utrans = {vb, Wb) in a least squares sense: 

'{-yi + Vh) 1 
(-2/2 + Vb) 1 

{xi - Xb) 

{X2 - Xb) 



OJ 
Vb 
Wb 



Vl 
V2 

Wl 
W2 



(17) 
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where b = {xb,yb)- 

With this we return to (5), calculating where ^ = ^ and ^ = Utrans in order to 

obtain a PDE similar to (6). We get 



db dTje)^ 

dt dt 



=r(6')u(x,t) + Utrans + 

Differentiating the rotation operator is straight-forward, 

d 



dt 



-X 



dt 



ne) 



-ft MO) 



CO 



Using this we get 



d 
dt 



cos(7r/2 + 6') -sin(7r/2 + 6l) 
sin(7r/2 + 6') cos(7r/2 + 6') 

= a;r(7r/2 + ^)- 



Mix) = r(^)u(x, t) + Utrans + CoT{Tr/2 + 0)x. 



(18) 



Inserting this into (5) gives a PDE, of the same type as (6), for the segment dynamics in the 
case when the velocity field decomposition (16): 



dl 
dt 



+ (r(0)Urest + Utrans + a;T(7r/2 + 6)^) 



dl 



1 



0. 



(19) 



The discretization of this PDE fits naturally into the explicit method summarized above, with 
two additions. First, apply Strang splitting (7) to obtain a second order accurate method. 
Secondly, the mapping parameters are updated before the coefficients {v,w) are reevaluated, 

gn+l = Qn_^ 
b"+l = b" + AtUtrans- 



3 Numerical results for interface tracking 

Here we provide numerical results for the pure interface tracking problem. So as to not 
introduce any new discretization errors, we let the convective velocity field be given by a 
closed expression (as opposed to computing an approximate solution to e.g. Navier Stokes 
equations) : 

U7(x, t) = cos(7rt) ( — cos(a;i — 7r/2) sin(a;2 — 7r/2) cos(x2 — 7r/2) sin(xi — 7r/2)) (20) 

U/7(x,t) = (sin(7ra:)^sin(27ry) — sin(7ry)^ sin(27rx)) (21) 

An example of a computation with u = u// can be seen in Figure 4. The curvature of the 
interface is known to to be difficult to compute, even to "eye norm" accuracy. It is even 
the case with many methods that some kind of filtering is needed to get acceptable curvature 
results, though this is not generally discussed. In Figure 4, we also show the interface curvature 
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at different times. The two sharp peaks are well captured, without any oscillations, and there 
are no visible artifacts at the segment boundaries. We compute curvature in local coordinates, 
via 

= (l + //(^)2)3/2' (22) 

as discretized with the usual centered difference approximations. 

For testing numerical convergence we construct test cases that have an analytical solution. 
The convective field u/ oscillates so that, at t = 1, the interface is returned to its initial 
configuration, which we take to be a circle. In Figure 5 wc present convergence results for 
both the interface location and its curvature as the grid is refined, with the explicit (LaxW) 
method. The time step is refined so that the ratio A^/At = 2 is kept constant. We see the 
expected second order convergence in both time and space for the interface, in ||x(r)— x(0) ||oo- 
In Figure 6 we give similar results for the implicit (CN) method. 

3.1 Convergence in curvature 

Also in Figure 5, to substantiate our claim that the interface curvature is well captured by 
our method, we give convergence results for curvature. Here we see first order convergence 
in \\k{s,T) — l/r||oo and close to second order convergence in \\k{s,T) — l/r||;2. That the 
2-norm convergence is substantially higher than in oo-norm tells us that the first order errors 
are highly localized. These results are with the LaxW method and we wish to emphasize 
that no filtering or smoothing of any kind were used to obtain them. Together with the 
second order convergence of the interface in oo-norm, we take these curvature convergence 
results as strong endorsement of the accuracy of the proposed method. In particular, the 
domain-decomposition approach seems vindicated. 

With the CN method we were not able to get equally convincing convergence results for 
the curvature. Still, in Figure 6, we have first order convergence of the curvature in 2-norm 
(in oo-norm we have convergence, but of order < 1). 

3.2 Mass, or area, conservation 

There is no inherent mass conservation in our method. On the other hand, there is no reason 
to expect it to suffer from poor mass conservation in the same way as level set methods are 
known to do. 

In Figures 5 and 6 (right column) we give mass loss results as the interface is refined, for 
the explicit and implicit method. That is, we compute the mass loss as a relative measure, via 
(mass(O) — mass(r))/mass(0) where the mass is simply computed as the area of the polygon 
enclosed by the interface. 

We note here that the mass loss shows second order convergence as the interface is refined, 
and that the constant is very small - several orders of magnitude below unity. 

3.3 Dynamic mapping parameters 

Finally, before we develop applications, we give numerical results for the extended method 
in Section 2.8, where the mapping parameters vary. Again with the oscillating circle case 
u = u/, we have second order convergence in the interface in oo-norm, first order convergence 
in curvature and second order convergence in the mass conservation, see Figure 7. In Figure 8 
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lO"* 10 ^ 10 ' 10 " 10 ^ 10 ^ 10 " 10" 

A^ 



Figure 5: Explicit (LaxW) method. Left: Convergence of interface in co-norm as grid is 
refined, with u = u/ (Eq. (20)), at t = 1. Middle: Convergence of interface curvature in (•) 
oo- and (*) 2-norm. Right: Convergence in relative mass loss (solid line denotes second order 
convergence) . 




A^ A^ A^ 



Figure 6: Implicit (CN) method. Left: Convergence of interface in oo-norm as grid is refined, 
with u = u/. Middle: Convergence of interface curvature in (*) 2-norm. Right: Convergence 
in relative mass loss (solid line denotes second order convergence). 
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Figure 7: Explicit (LaxW) method, with velocity field decomposition (Section 2.8). Left: 
Convergence of interface in oo-norm as grid is refined, with u = u/. Middle: Convergence of 
interface cm'vature in (•) oo- and (*) 2-norm. Right: Convergence in relative mass loss (solid 
line denotes second order convergence). 




-1 -0.5 0.5 1 0.4 0.45 0.5 



Figure 8: Left: Interface at four times, and pure rotation velocity field. Left: Close up of 
boxed region in left figure, with interface at t = as solid red line, and interface at i = 27r as 
dots. These are indistinguishable. 

we give results where the interface is advected in a pure rotation velocity field, as is a common 
test case for interface tracking methods. A notable characteristic of our method is that, aside 
from initialization errors, there are no numerical errors in this computation, because of the 
velocity field decomposition. 

4 Boundary integral method for Stokes in 2D 

When applicable, boundary integral methods are often efficient. Here we give a boundary 
integral formulation for Stokes flow in 2D using our method. This illustrates some important 
simplifications that arise when using our method. 

For an introduction to the field of boundary integral methods for viscous flow, we refer to 
the textbooks Pozrikids [32] and Kim and Karrila [33], as well as a more recent survey paper 
by Pozrikidis [34]. 

Using the so called single layer formulation, the velocity at any point, xq G ft C M", can 
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be evaluated as 

Uiixo) = j^Gij{x,xo)qj{x)dS{x), (23) 

where Gij is a Greens function (fundamental solution) called the Stokeslet, and qj is a source 
density on the interface. Here we consider a free space problem in M^. Then the Stokeslet is 
given by 

Gij{x,xo) = -6ijlog{r) + 

where X = X — xq, r = \x\ and Sij is the Kronecker symbol. There remains only to compute 
the densities qj, which in general requires solving a (full) linear system. However, in the case 
where the viscosity inside is equal to the viscosity outside of the interface, one obtains 

where s parametrizes F and is the viscosity. These results were given in [35]. The force 
acting on the interface is simply 

fj{s) = aK{s)hj{s), 

where a is the coefficient of surface tension, k is the curvature and hj are the unit normals. 
In terms on non-dimensional variables, we write 

Qj{s) = 

where Ca := ^ is known as the capillary number (for some typical velocity scale U). 

Now we pose the evaluation of (23) in the setting of our segments. By definition we have 
that any integral over the interface can be split into a sum of integrals over the segments 7^, 

/ g{s)ds = Y.[ mm (24) 

In order to perform this change of variables in (23) from a global parametrization, s, to local 
variables (C, /('?)) we note a few invariant properties of the integrand. First, the curvature is 
a geometric invariant, so instead of computing k{s) we can compute it in local coordinates 
directly, as given by (22). Similarly, for the normal vectors we directly compute 

n{0 = [f'{0 -!]• 

Finally, note that the global-to-local mapping, M. (cf. Eq. (3)) is distance preserving. Hence, 
we directly express the Stokeslet in local coordinates 
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where x := (C, /(C)) — Mxq and r = ||x||. For a point xq away from the interface, using the 
inverse mapping (4), we have that 

^(^o) = -^l]A_/^^«(OG'.,(C,/(0,^o)n,(C)dC/ 

= -i^E a/^^ ^^G.(C,/(e),xo)n,(OdC (25) 

1^ k '>''= 

These integrals arc naturally discretizcd on the segment grids, where the start of the inte- 
gration interval, denoted C', corresponds to the end of the previous segment, i.e we compute 
integrals of the form 

These are discretized with trapezoidal quadrature over the ^-grid. 



4.1 Treatment of singular integrands 

As the point xq approaches the interface the integrands in (25) become singular. In order to 
evaluate Uj on the interface (e.g. for time-dependent problems) we need to provide precisely 
how to compute these integrals - showing that they exist and what the relevant limits are. 

Suppose xo lies on the interface. On some segment, xq corresponds to a point in local 
coordinates (Co, /(Co))- The distance function then takes the form 



r = r(C, Co, /) = V(C-Co)2 + (/(C)-/(Co))2. 



(26) 



First, we have the log-type singularity due to the first term in G, which needs careful 
treatment. Secondly, the terms involving XiXj/r^ are non-singular, but the limits as C — >■ Co 
need to be determined. 

The logarithmic singularity can be subtracted away, in two steps: 



Uj 



g,(C)log(r(C,Co,/))dC 
f-6 



= / fe(C)-g,(Co)]log(r(C,Co,/))ciC + «,(Co) / log(r(C,Co,/))dC 

J a J a 

= rfc(C) -Q,(Co)]log(r(C,Co,/))dC+ 

Ja 



+ 9i(Co) 



(27) 



Now we show that the first two integrands are non-singular, compute the corresponding limits 
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as ^ — )• ^0 and compute the last integral analytically. For the latter we get 



J a 



(28) 



(log(|6- xol) - xo) ifa = a;o 

-(log(|a - xo|) - l)(a - xo) if 6 = xq 

(log(|6 - xo|) -l){b- xo) - (log(|a - xo|) - l)(a - xq) otherwise. 

The integrand in the first integral in (27) behaves as xlogx for small x, and hence goes to 
zero. More specifically, the limits of the two integrands (i.e. for the x- and y-components of 

u) are 



lirn {q({0 - q(m log(r(6) = Hm 



and 



lirn (g|(0 - ^^'(eo)) log(r(0) = hm 



J I 



-no 



«^«o Vl + /'(0' i + /'(eo) 



logr(0 = 0. 



For the second integral in (27) the relevant limit is 



l + f'H^o)>0, 



so that the second integrand goes to log(\/l + /'^(^o)) near ^ = ^o- This concludes the 
treatment of the integrals involving log-singular integrands. The simple quadrature rule 
above is applied and the limit values are used when the integrands are not immediately well 
defined. 

Similarly, for the terms in G involving XjXj/r^ we provide the required limits. 



lim 



lim 



^^^0 r{0^ l + /'(Co)2 

(e-Co)(/(0-/(Co)) nco) 



lim 



im) - mo))' 



/'(eo)^ 



i + /'(eo)2- 

This concludes the formulation of the boundary integral method. 



4.2 Numerical results for BIE Stokes 

The boundary integral method above fits into the proposed interface tracking method, simply 
by letting (25) define the convective field, u, in e.g. the explicit (LaxW) time step scheme. 

In Figure 9 we give illustrative plots of the flow field, evaluated over a grid, as computed 
via the boundary integral (25). These show the expected flow features for a free space flow 
driven by the curvature of the interface. For a comparison to Navier-Stokes flow, and results 
for Stokes for a shear flow, see Section 5.5.3 and Figure 11. Note that evaluating the velocity 
field away from the interface is only done for post processing, e.g. visualization. The resolution 
here was = 1/100, and the interface was initialized with four segments (as in Figure 1). 



19 




Figure 9: Interface in free space Stokes flow, driven by surface tension, at four times. The 
convective field u is obtained by evaluation the (singular) integrals in (25) on the interface. 
The flow fields shown are evaluations of (25) on a grid, used for visualization only. 
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10"® 10"" 10"^ 10"^ 10"^ 

Figure 10: Convergence of boundary integral method (25), for u and v, with second order 
given by sohd Hne. Left: random sample set of observation points xq away from the interface. 
Right: Observation points on the interface, i.e. singular integrands. 



The quadrature method described is trapezoidal, but with one interval of unequal length. 
In Figure 10 we establish second order convergence in oo-norm away from the interface and 
second order convergence in /i-norm on the interface, as the interface is refined. Here we have 
taken the interface circular, so that, by symmetry, the exact solution is u = Vxq. Away 
from the interface we randomize a set of observation points, and on the interface we take as 
observation points the entire interface. 



5 Incompressible Navier Stokes flow computations 

We now demonstrate that our interface tracking method is applicable to the well established 
application of two-phase Navier-Stokes flow and that good results are obtained. There is not 
much here which is specific to our method. Indeed, this is an important point - our interface 
tracking method works well with established components in multi-phase flow computations. 

5.1 Overview 

The task of solving a multi-phase Navier-Stokes flow system is complicated for several reasons. 
At a basic level, solving the Navier-Stokes equations themselves is a challenging problem on 
its own. The incompressible Navier Stokes equations are given by 

p(^^ + (u-V)u^ =-Vp + /iV2u + F (29) 

V • u = (30) 
-|-(Boundary conditions). (31) 

where denotes viscosity and p density. We will, for the remainder of this paper, consider 
the momentum equation in non-dimensional form: 

^ + (u-V)u = -Vp+— V^u + F, (32) 

with Re := the Reynolds number. 
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The two main difficulties here are that the first (i.e. momentum) equation is non-Hncar, 
and that the space of admissible solutions is restricted (by the second equation) to be diver- 
gence free in the entire domain. We refer the reader to the excellent textbook by Acheson 
[36]. Several classes of numerical methods exist for Navier-Stokes, based on finite volume 
(FV), finite difference (FD) and finite element (FEM), as well as spectral methods. These 
often have high accuracy orders with respect to the spatial discretization. There are, however, 
much fewer methods that are fully second order accurate also in time. In the fiuid mechanics 
community, this is generally regarded as of little consequence, since there are often physical 
time scales in the fiow that are small. In this work we use a fairly basic pressure correction 
FD method, as presented in Section 5.2. 

When a dynamic interface is added, e.g. the case of two immiscible liquids, the picture 
becomes more complicated still. The physical parameters of the fiow may be different in 
the different fiuids, and a two way coupling is introduced between the interface and the fiow 
field. That is, the fluid acts upon the free surface which simultaneously acts upon the fluid. 
From a mathematical point of view it is clear that the equations of the two (or more) fluids 
and the dynamics of the free surface constitute one large, coupled and highly non-linear 
system. Furthermore, application specific challenges arise frequently. In the ubiquitous case 
of a bubble in a bulk fluid the most widely used formulation introduces a singular source 
term in the momentum equation (29) to mediate forces on the fluid from the interface, e.g. 
surface tension. The lack of regularity in the velocity and pressure solutions to the Navier 
Stokes equations presents significant numerical challenges. Here, we shall use the immersed 
boundary (IBM) and immersed interface (IIM) methods as presented in the following sections. 



5.2 Numerical method for Navier Stokes 

The numerical method to treat the Navier Stokes equations (29)-(31) is a finite difference 
method on a staggered (MAC) grid, which employs so called Chorin splitting. It is presented 
in detail in the textbook by Strang [37, Ch. 6.7], and a basic implementation due to Seibold 
is available on-line [38]. We summarize the time-discrete method here: 

Compute a first intermediate velocity field, u*, by treating the non-linear term explicitly, 

"^At"'' ^ (u"- V)u" + F. (33) 

To avoid a severe time-step restriction the diffusive term is treated implicitly and a second 
intermediate is obtained via 

** _ * 

, = ^^Au". (34) 
At Re ^ ' 

This, together with a discrete approximation for the spatial derivatives, constitute a linear 
system. Several appropriate solvers exist for this system depending on the memory constraints 
and the boundary conditions applied. The final result is obtained by applying a correction 
step. 



u 



n+l 



= -Vp"+\ (35) 



At 

based on the pressure, which is obtained by solving a Poisson equation. 



-Ap"+^ = -^V • u**. (36) 
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We use central difference approximations for the spatial derivatives, so that V is replaced 
by the standard discrete nabla operator on a staggered grid Vh- However, when applied to 
the convection step a moderate upwind bias is applied so as to not introduce a numerical 
instability (see [37, Ch. 6.7]). 

Several remarks are in order. The pressure correction strategy is by no means the only 
one available for the task of decoupling the momentum equation from the incompressibility 
constraint. Another popular approach is the ^-splitting method (see e.g. Glowinski and 
Pironneau [39]). A large body of work exist on the topic of pressure correction methods 
for Navier Stokes. The approach here dates back to Chorin [40] and may be seen as the 
most elementary method. Several more recent methods were presented and evaluated in a 
recent survey by Guermond et. al. [41]. The class of methods that Guermond refers to as 
consistent splitting methods can be seen to be more accurate than Chorin sphtting. We also 
remind the reader that pressure correction methods do not necessarily have the property that 
both the non-divergence constraint and the prescribed boundary conditions are accurately 
satisfied simultaneously (due, at least in part, to the lack of a consistent numerical boundary 
condition for the pressure). Shedding light on these interesting topics is beyond the scope of 
this work, and we also content ourselves with using a first order accurate time discretization. 
A final remark concerns the efhciency of the method above: the linear solves in the diffusion 
and pressure correction steps involve constant system matrices. These matrices need only be 
factorized once, leaving just the back- and forward substitution steps to be performed in the 
time loop. The major constraint here is clearly the memory requirements of the factorized 
matrices, which may be significantly less sparse, but the efficiency obtained is of great practical 
value. 

The flow solver couples to the interface tracking method previously presented, using Strang 
splitting (7). 

5.3 Interface/flow coupling: immersed boundary formulation 

As a basic case, let F be a closed interface separating two liquids with equal density and 
viscosity. Two main approaches exist for coupling the motion of a dynamic interface with a 
incompressible, viscous bulk flow (described by the Navier Stokes equations above). First, 
separated grids may be generated for the separated flow regions - so called body fitted meth- 
ods. The second class of methods follows the work of Peskin and collaborators [42, 43, 44, 45] 
on the immersed boundary (IB) method, where one views the separated flow regions as one 
system and a suitable force density is added to mediate the action of the interface on the fluid 
system. This method dates back to the 1970's but still sees a remarkable amount of use. The 
2002 survey paper by Peskin may be of particular interest to some readers [45] . 

Here, we shall use the IB formulation, with surface tension driving the flow internally. 
It is well known that the surface tension force as a function of Lagrangian variables on the 
interface can be formulated as 

{{s,t) = aK{s,t)h{s,t), (37) 

where a G is the surface tension coefficient, k is the interface curvature and n is the unit 
inward normal. Recall that the capillary number, Ca, is defined as Ca := We may then 
write the non-dimensionalized force that will enter the momentum equation (32) in Eulerian 
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variables as 



F(x, t) = I Kis, t)his, t)Srds, (38) 

where is a Dirac measure on the interface, i.e. (5r(x) = 5(x — r(s)) and 5(x) is the 

composition of one-dimensional delta functions, 5(x.) = S{xi)S(x2). Transformations from 
Lagrangian to Eulerian variables in this way are at the heart of the IB method and much 
related work. It is also appropriate to precisely define how the convective field in the interface 
advection ODE (1) arises: 

u(r,i) := / u(x,i)Mx. (39) 
Jn 

In order to complete the formulation of the discrete immersed boundary method one 
needs to introduce discrete delta function approximations. We refer the reader to a study by 
Tornberg and Engquist on the treatment of singular source terms [46] . We use the piecewise 
cubic delta function approximation introduced therein, 

M.j^i*/*'- (40, 
1^0, \x\ > 2h 

with 

, fl-iM-lrf + kM'', 0<r<l 

This delta function has narrow support on the grid and possesses a larger number of discrete 
moments than other commonly used variants. There are other good choices for the discrete 
delta function approximation. However, there is a generally held view, concerning stability, 
that the same approximation should be used in the "spreading" step (38) as in the evaluation 
step (39). 



5.4 Immersed interface method 

The singular source term introduced by the IBM gives rise to a jump discontinuity in the 
pressure, which is well understood physically: the pressure gradient balances the interfacial 
forces. This jump is 

[p]{s,t)=f{s,t)-h{s,t). (42) 

That is, the pressure jump equals the magnitude of the interfacial force in the normal direc- 
tion, which is simply [p](s, i) = aK{s,t) in our case. Additional jump conditions need to be 
considered in settings where there are e.g. tangential stresses on the interface. 

The immersed interface method (IIM), introduced and developed by LcVcque and collab- 
orators in a series of papers [47, 48, 8], is in essence a evolution of the IBM where the jump 
conditions are corrected for in the numerical method rather than taken into the momentum 
equation as a source term. The textbook by Li and Ito [9], and a recent paper by Le et. al 
[10] may also be of interest. Derivations of the correction terms are given in the works cited 
above, and we humbly restate the relevant results here. 
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Going back to the NS method (33)-(36), drop the source term, F completely from the 
convective step, (33). The diffusive step (34) is unaltered, and the pressure correction step 
(35) and (36) we replace by 




-Ap"+i = -^V-u** + C. 



(43) 



The correction terms, B = {B^, B^) and C, only have meaningful definitions as discrete 
functions. Recall that we are working on a staggered (MAC) grid. It turns out to be natural to 



This concludes the formulation of the discrete immersed interface method. 

5.4.1 IBM vs. IIM: Accuracy 

Due to the smearing of the interface introduced by the regularized delta function, the IBM 
will at most attain first order accuracy for problems with non-smooth solutions. The IIM has 
been shown to produce second order accurate results. 

5.4.2 IBM vs. IIM: Stability 

Interestingly, there has been much recent work on the stability of the IBM, despite the many 
years since it was introduced. Peskin himself emphasizes stability as one of the outstanding 
challenges for the IBM in his Acta Numerica paper [45]. Newren et. al. [49] appear to 
have resolved this issue by providing a general stability theory for the IBM applied to Stokes 
flow. They show unconditional stability for a Crank-Nicholson method and backward Euler 
method. Following that paper, Newren et. al. [50] discuss practical solver strategies for 
semi- and fully implicit immersed boundary methods for Navier Stokes. They substantiate 
the widely held view that fully explicit methods are competitive, despite the short time steps, 
due to the high computational costs associated with the iterative methods needed in implicit 
IBM formulations. In another recent paper, Hou and Shi [51] also obtain unconditionally 
stable discretizations for the Stokes case. Based on earlier work by Hou, Lowengrub and 
Shelley [52] and their unconditionally stable method, they go to great lengths to derive an 
efficient semi- implicit method with good stability properties. While this appears to be the 
most efficient method known at present, it relies heavily on spectral properties of the spatial 
discretization and, hence, is only applicable in certain cases. 

Newren et. al. point out two conditions as sufficient to prove stability of immersed 
boundary methods: First, the spreading and interpolation must be adjoint operators. Second, 
the velocity field needs to be discretely divergence free. Hou and Shi reiterate this conclusion. 





(44) 



(45) 
(46) 



Cij = (V/i • B)y. 
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Hou and Shi also make a separate point about the need to aheviate the stiffness in the 
interface/flow couphng, and go about it by using the arc-length and tangent angle formulation 



Little is known to date about the stability properties of immersed interface methods. There 
are several vague comments on stability in the IIM literature [48, 8, 10, 9]. The authors seem 
to agree that the explicit methods suffer from a severe stiffness, which motivates implicit 
methods. Going back to the sufficient conditions for stability of IB methods, it is clear 
that there is no adjointness property applicable in the IIM case (since there is no spreading 
operator). However, one could reasonably expect the non-divergence condition to be better 
handled by the IIM: In the IBM, as the pressure jump gets bigger the smeared pressure 
becomes under-resolved numerically and oscillatory near the interface. This would lead to a 
less accurate pressure correction, and hence a failure to be accurately divergence free in the 
velocity. The pressure solve is more accurate in the IIM. Further, the idea of formulating the 
interface problem in terms of arc-length and tangent angle has not been applied to the IIM 
as far as we know. 

5.5 Numerical results: Navier-Stokes 
5.5.1 Shear flow 

We illustrate the NS/IBM solver by considering a drop in shear flow. The shear velocity is 
imposed as a Dirichlet boundary condition on the top and bottom boundaries, ramped up via 



for some ramp-up time too and shear velocity -Uoo- We let ^(0, t) = —u{l,t). In the x-direction 
we impose periodic boundary conditions. In Figure 11 wc present computational results with 
different coefficients of surface tension. Here WC used = 1, too = 1/5 and 80 grid-points 
in the y-direction. These results could equally well have been obtained with the immersed 
interface method given above. Throughout this section we use the channel height, Ly = 1, 
when determining the Reynolds number. 

5.5.2 Convergence 

The method used here is not more than first order accurate, as discussed prior. To verify this 
we ran numerical convergence tests using a shear flow setup. Here we use an initially circular 
drop with radius r = 1/3 centered in a unit square domain. The shear boundary conditions 
arc ramped up as above, with too = 1/5 and we compute until t = 1/4. With Re = Ca = 1 we 
get, at that time, that the interface is moderately deformed, but not close to a steady state. 

In Tables 1 and 2 we give convergence results for the immersed boundary and immersed 
interface methods respectively, as the grid is refined. The ratio h/At is kept fixed, so several 
hundred time steps are taken on the finest mesh. These results indicate the expected overall 
first order convergence for both methods. 

Here we let Ui be a sequence of computed solutions, and assume the error with respect to 
some unknown exact solution u behaves like \\ui — u\\ ^ 0{h^) + O(AtP). We define A* := 
\\ui — Uj-ill for various norms and have estimates of the convergence order p = log2(A*"'"-^/A*) 
as the grid is refined by a factor 2. To see a benefit in terms of accuracy from using the 



(cf. [52, 51]). 
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Figure 11: Interface location at various times, and flow field at t = 8. Top pair: Re = 1, Ca = 
1. Below critical Ca, the interface reaches a steady state. Bottom pair: Re = l,Ca = 2. 
Above critical Ca, no steady interface is obtained. 
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Table 1: NS/IBM grid refinement study, shear flow, Top: li, Middle: I2, Bottom: l^. See 
Section 5.5.2 

immersed interface method, a second order accurate time stepping method for the flow solver 
would have to be implemented. 

5.5.3 Brief comparative study 

Whereas a full comparison between the immersed interface- and boundary method including 
solution methods would be valuable, it is beyond the scope of this work. The space of 
parameters, e.g. discretization and solution methods, capillary- and Reynolds numbers, is 
very large. Still, wc think some brief comparisons arc in order. 

First, for the simple shear flow setup described above we compare the results obtained 
with the immersed interface- and boundary methods with explicit time stepping. As can be 
seen in Figure 12, there is a small but visible discrepancy. Here, the domain was [0 2] x [0 1], 
with a grid with hx = hy = 1/100 and the interface was discretized with = 1/200. Again, 
the physical parameters were Re = Ca = 1 . In this test we noted that the IIM was more prone 
to going unstable, as expected. We had to take At = 1/300 with the IBM and At = 1/3000 
with the IIM, to get stable computations. 

For a more challenging case, e.g. with Ca several orders below unity, we expect to see the 
implicit (CN) method stable for larger time steps than the explicit method. However, this 
has not been evident to us - at least not to the extent that the overall computational cost of 
a run becomes smaller with an implicit method. That implicit methods fail to be competitive 
for the IBM is in line with results by Newren et. al. [50] and others. We have to conclude 
with noting that more research is need to illuminate the eflicient deployment of immersed 
boundary- and interface methods. 

In Figure 13 we consider a shear flow case where the Reynolds number goes to zero, i.e. the 
limit of Stokes flow. Here we compare the results from the Navier-Stokes (IBM) method with 
the boundary integral Stokes method (Section 4) and find excellent agreement as Re ^ 0. 
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Table 2: NS/IIM grid refinement study, shear flow, Top: li, Middle: I2, Bottom: /qo- See 
Section 5.5.2 
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Figure 13: Interface in shear flow, at time T = 1, and Ca = 5. We illustrate that as Re — >■ 0, 
the Navier-Stokes (IBM) method converges to the results from the boundary integral Stokes 
method (cf. Section 4). Here the interface is initially circular, with diameter 1, and at location 
(3,3/2), in a domain [0,6] x [0,3] and shear boundary conditions discussed in Section 5.5.1. 
Note that the T = 1 is before a steady state is reached. Bottom: close-up from above. 

6 Conclusions and outlook 

We have introduced an interface tracking method which is both Eulerian and explicit, as 
motivated by current and future needs for micro- and complex flow modeling and simulation. 
For instance, the PDE-based methods by Khatri and Tornberg [6] for simulating soluble and 
insoluble surfactants fit naturally into this framework. 

The method presented is based on a domain decomposition approach, where the interface 
is split into segments, or patches. On each patch a hyperbolic PDE is solved in time to update 
the interface location, and adjacent segments are coupled via an exchange of ghost points. 
One explicit and one implicit time step method was presented, both second order accurate. 

Numerically, the accuracy and consistency order of the proposed method was established: 
second order in oo-norm for the interface location. We were also able to present first and 
second order convergence in the interface curvature - a quantity which is known to be hard 
to compute (without filtering). Mass conservation was found to be well behaved, converging 
at second order with a constant several orders of magnitude below unity. 

An extension of the method was also presented, which made use of a decomposition of 
the velocity field into translation and rotation components, as motivated by a variety of flow 
applications where such components tend to dominate. The resulting method was also shown 
numerically to have the desirable accuracy characteristics mentioned above. 

We treated two applications: Stokes and Navier-Stokes flow. The first used a boundary 
integral method, and illustrated how an application with integration over the interface can 
be treated in our framework. We found that it was simple to derive and analyze a method 
which was second order accurate for integrals of the Greens' function for Stokes flow in free 
space. 
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Finally, with the application of our method to two-phase incompressible Navier Stokes 
flow, we showed that our method works well with existing components for such problems. 
Both the immersed interface- and immersed boundary method was used. We gave and eval- 
uated a first order accurate method for solving the Navier-Stokes equations and presented 
computations of a bubble in shear flow. The two-phase flow results are by no means new, but 
serve to underline that our method can be used in established applications with good results. 

However, as of yet, we have not generalized our method to 3D in sufficient generality. The 
framework and equations generalize naturally - the difficulty is representing closed surfaces 
as a union of patches, in a way that allows the patches to couple together robustly. One would 
like to use regular grids in 2D for these patches, to be able to use well-established numerical 
method for advecting the interface, computing curvature and quadrature etc. 

None the less, we believe that at a method that is explicit in the interface location and 
Eulerian is very useful for treating complex interface tracking applications. There is a wealth 
of interesting flow applications on the micro scale where continuum modeling on the interface 
becomes relevant. 
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